!/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!*  Copyleft (!c) 2004-2005 Div. Nucl. Med., Cyric, Tohoku Univ.    *
!*  Copyleft (!c) 2005-2011 Elseif Laboratory.                      *
!*  elseifkk@users.sf.net                                           *
!*                                                                  *
!*  All Rights Reversed.                                            *
!*                                                                  *
!*                                                                  *
!*  LICENSE AGREEMENT:                                              *
!*                                                                  *
!*  You agree that we don't agree to anything.                      *
!*                                                                  *
!*                                                                  *
!*  DISCLAIMER OF WARRANTIES:                                       *
!*                                                                  *
!*  This software is provided to you "AS IS," and the authors and   *
!*  contributors disclaim any warranty or liability obligations to  *
!*  you of any kind, implied or express, including but not limited  *
!*  to the warranties of merchantability, fitness for a particular  *
!*  purpose, and non-infringement.                                  *
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
module Ground_Control
  use muty
  use cuty
  use buty
  use futy
  use kuty
  use optimizer
  use MatchMaker
  use ana
  implicit none
  !-----------------!
  integer*4 mana    ! 
  !-----------------!
  character*(*),parameter::BUILD_TIME="Build on "//"@BUILD_TIME"//"\n"
  character*(*),parameter::PROGRAM="\nana version 0\n"//&
       "Copyleft (!c) 2004-2005 Div. Nucl. Med., Cyric, Tohoku Univ.\n"//&
       "faely version 0.0.1\n"//&
       "Copyleft (!c) 2005-2011 Elseif Laboratory."
  character*(*),parameter::HELP=&
       "Usage:\n"//&
       "faely [Options] target [reference] [output]\n\n"//&
       "Options:\n"// &
       "  -x<number>,-y<number>,-z<number>\n"// &
       "    set translation vector.\n"//&
       "  -1<number>,-2<number>,-3<number>\n"// &
       "    set Euler angle.\n"//&
       "  -n<number>,-m<number>,-k<number>\n"// &
       "    set reduction factor.\n"// &
       "  -R<number>\n"// &
       "    set reduction factor for all.\n"//&
       "  -X<number>,-Y<number>,-Z<number>\n"// &
       "    set voxel size.\n"//&
       "  -L<number>\n"// &
       "    set cubic voxel size.\n"//&
       "  -s<number>\n"// &
       "    set median filter size in x-y plane.\n"//&
       "  -S<number>\n"// &
       "    set median filter size.\n"//&
       "  -l<number>,-u<number>\n"//&
       "    set lower and upper threshold.\n"//&
       "  -C<number>\n"//&
       "    set number of cluster.\n"//&
       "  -D<number>\n"//&
       "    set output datatype; available=2,4,8,32,64.\n"//&
       "  -d<number>\n"//&
       "    set internal datatype; available=2,4,8,32,64.\n"//&
       "  -i\n"//&
       "    evaluate mutual info between input and output.\n"//&
       "  -P\n"//&
       "    evaluated information-potentialed image and all.\n"//&
       "  -e\n"//&
       "    use histogram equalization.\n"//&
       "  -H\n"//&
       "    set number of hist bin.\n"//&
       "  -r<number>\n"//&
       "    linear remapping of histgram from 0-<number>.\n"//&
       "  -p<number>\n"//&
       "    set pdf smoothing width.\n"//&
       "  -b<number>, -B<number>\n"//&
       "    use adaptive binary filter.\n"//&
       "  -G<number>\n"//&
       "    cut background.\n"//&
       "  -I<interpolation method={LIN|NN|SINC|PV}>\n"//&
       "    chose interpolation method for internal calculation.\n"//&
       "  -M<interpolation method={LIN|NN|SINC}>\n"//&
       "    chose interpolation method for creating transformed image.\n"//&
       "  -c\n"//&
       "    exchange tar->ref.\n"//&
       "  -t\n"//&
       "    save moved image.\n"//&
       "  -T\n"//&
       "    save moved image in the same dimension as the reference.\n"//&
       "  -Q\n"//&
       "    do not make output.\n"//&
       "  -f\n"//&
       "    ignore calibration factor.\n"//&
       "  -+\n"//&
       "    ignore negative when truncated-read.\n"//&
       "  -E\n"//&
       "    Extract connected components.\n"//&
       "  -g<number>\n"//&
       "    Gaussian smoothing(3D) width in mm\n"//&
       "  -o\n"//&
       "    set output filename.\n"//&
       "  -W\n"//&
       "    search R and D at the same time.\n"//&
       "  -a\n"//&
       "    allow elastic transform for registration.\n"//&
       "  -F<number>\n"//&
       "    selects input file format. 1=shimadzu.\n"//&
       "  -U\n"//&
       "    multiply tar and ref.\n"//&
       "  -w\n"//&
       "    for experiment"
  integer*4,parameter::NOP                   = Z'00000000'
  integer*4,parameter::VERBOSE               = Z'00000001'
  integer*4,parameter::REDUCE_R              = Z'00000002'
  integer*4,parameter::REDUCE_V              = Z'00000004'
  integer*4,parameter::REDUCE_AO             = Z'00000008' 
  integer*4,parameter::REDUCE                = Z'0000000E' ! mask
  integer*4,parameter::MOVE                  = Z'00000010'
  integer*4,parameter::USE_MEDIAN_FILTER     = Z'00000020'
  integer*4,parameter::USE_MEDIAN_FILTER_XY  = Z'00000040'
  integer*4,parameter::MASK                  = Z'00000080'
  integer*4,parameter::CLUSTERIZE            = Z'00000100'
  integer*4,parameter::GET_MI                = Z'00000200'
  integer*4,parameter::GET_IP                = Z'00000800'
  integer*4,parameter::USE_HISTEQ            = Z'00001000'
  integer*4,parameter::USE_REMAP             = Z'00002000'
  integer*4,parameter::DO_EXP                = Z'00004000'
  integer*4,parameter::REF_READ              = Z'00008000'
  integer*4,parameter::USE_ADAPBIN_FILTER    = Z'00010000'
  integer*4,parameter::USE_ADAPBIN_FILTER_XY = Z'00020000'
  integer*4,parameter::CUT_BG                = Z'00040000'
  integer*4,parameter::DUMP_PDF              = Z'00080000'
  integer*4,parameter::SAVE_MOVED            = Z'00100000'
  integer*4,parameter::EXCHANGE_INPUT        = Z'00200000'
  integer*4,parameter::NO_OUTPUT             = Z'00400000'
  integer*4,parameter::EXTRACT_CCS           = Z'00800000'
  integer*4,parameter::USE_GSMT              = Z'10000000'
  integer*4,parameter::ELASTIC               = Z'20000000'
  integer*4,parameter::COUPLED               = Z'40000000'
  integer*4,parameter::MULTIPLY              = Z'80000000'
  integer*4,parameter::ADF_MASK              = Z'03000000'
  integer*4,parameter::IGNORE_CF             = Z'01000000'
  integer*4,parameter::IGNORE_NEG            = Z'02000000'
  integer*4,parameter::CHG_EXTENSION         = Z'04000000'
  integer*4,parameter::MATCH_DIMENSION       = Z'08000000'
  !
  integer*4,parameter::MANA_MUL_IMG = 1
  !
  integer*4,parameter::PATH_LEN_MAX = 256
  integer*4,parameter::FF_ANALYZE = 0
  integer*4,parameter::FF_SHIMADZU = 1
  type cmd_info
     real*8 dx(3)
     real*8 dt(3)
     real*8 dm(3)          !! internal use 
     real*8 ds(3)          !! internal use
     real*8 red(3)         !! reduction factor
     real*8 d(3)           !! voxel size for reduction
     real*8 gsmt_s         !! gaussian smoothing's sigma in mm
     real*8 bg_level       !! background reduction factor
     real*8 lower,upper    !! threshold
     real*8 pdf_smtw       !! width of gaussian smoothing pdf
     integer*4 nc          !! number of cluster
     integer*4 n_mf(3)     !! median filter size
     integer*4 n_abf(3)    !! adapbin filter size
     integer*4 int_kind    !! data type internal
     integer*4 out_kind    !! output data type
     integer*4 int_type_in !! interpolation method
     integer*4 int_type_mv !! output interpolation method
     integer*4 nh          !! number of hist bin.
     integer*4 ns          !! number of gray scale
     integer*4 ff          !! input file format
     character(PATH_LEN_MAX) fin  ! not to be changed
     character(PATH_LEN_MAX) fref
     character(PATH_LEN_MAX) fout
  end type cmd_info
  private print_help
  private print_version
contains 
  !------------------------------------------------------!
  subroutine print_help
    implicit none
    call messn(HELP)
  end subroutine print_help
  !------------------------------------------------------!
  subroutine print_version
    implicit none
    call messn(PROGRAM)
    call messn(BUILD_TIME)
  end subroutine print_version
  !------------------------------------------------------!
  subroutine init_ci(ci)
    type(cmd_info),intent(out)::ci
    ci%dx=0
    ci%dt=0
    ci%dm=1
    ci%ds=0
    ci%fin=""
    ci%fref=""
    ci%fout=""
    ci%red=1
    ci%d=0
    ci%pdf_smtw=0
    ci%n_mf=0
    ci%n_abf=0
    ci%bg_level=0
    ci%nc=0
    ci%int_kind=DT_RKIND_8  ! real8
    ci%out_kind=DT_UNKNOWN  ! same as input
    ci%lower=-HUGE(0.d0)
    ci%upper=HUGE(0.d0)
    ci%nh=0
    ci%ns=0
    ci%int_type_in=0
    ci%int_type_mv=0
    ci%gsmt_s=0
    ci%ff=0 ! analyze format
  end subroutine init_ci
  !------------------------------------------------------!
  subroutine me_first_mana(ci)
    integer*4 i,j,n,len
    character*128 argv
    real*8,pointer::x
    integer*4,pointer::k
    character(PATH_LEN_MAX),pointer:: s
    type(cmd_info),intent(out),target::ci
    integer*4 istat,cd
    character*32 cnt
    logical*4 do_cycle
    integer*4,pointer::int_type
    n=command_argument_count()
    mana = NOP
    istat=0
    call init_ci(ci)
    cd=0
    if(n==0) then
       call print_version
       call print_help
       call exit(0)
    end if
    nullify(x,k,s)
    do i=1,n
       if(green(mana, VERBOSE)) then
          if(cd==0) then
             call mess("<<< Ground control to mana, "//&
                  "commencing count down, engines on..., ")
             cd=1
          else
             call mess(cnt(1:d2str(n-i+1,cnt))//"..., ")
          end if
       endif
       call get_command_argument(i,argv,len)
       if(read_arg()) cycle
       if(argv(1:1)=='-' .and. len>1) then
          select case(argv(2:2))
          case('F')
             k=>ci%ff
          case('g')
             mana = ior(mana, USE_GSMT)
             x=>ci%gsmt_s
          case('G')
             mana = ior(mana, CUT_BG)
             x=>ci%bg_level
          case('C')
             k=>ci%nc
             mana = ior(mana, CLUSTERIZE)
          case('H')
             k=>ci%nh
          case('d')
             k=>ci%int_kind
          case('D')
             k=>ci%out_kind
          case('l')
             x=>ci%lower
             mana = ior(mana, MASK)
          case('u')
             x=>ci%upper
             mana = ior(mana, MASK)
          case('x')
             x=>ci%dx(1)
             mana = ior(mana, MOVE)
          case('y')
             x=>ci%dx(2)
             mana = ior(mana, MOVE)
          case('P')
             mana = ior(mana, GET_IP)
          case('p')
             x=>ci%pdf_smtw
          case('z')
             x=>ci%dx(3)
             mana = ior(mana, MOVE)
          case('X')
             x=>ci%d(1)
             mana = ior(mana, REDUCE_V)
          case('Y')
             x=>ci%d(2)
             mana = ior(mana, REDUCE_V)
          case('Z')
             x=>ci%d(3)
             mana = ior(mana, REDUCE_V)
          case('L')
             x=>ci%d(1)
             mana = ior(mana, REDUCE_V)
          case('1')
             x=>ci%dt(2)
             mana = ior(mana, MOVE)
          case('2')
             x=>ci%dt(1)
             mana = ior(mana, MOVE)
          case('3')
             x=>ci%dt(3)
             mana = ior(mana, MOVE)
          case('n')
             x=>ci%red(1)
             mana = ior(mana, REDUCE_R)
          case('m')
             x=>ci%red(2)
             mana = ior(mana, REDUCE_R)
          case('r')
             k=>ci%ns
             mana = ior(mana, USE_REMAP)
          case('k')
             x=>ci%red(3)
          case('I','M')
             if(argv(2:2)=='I')then
                int_type=>ci%int_type_in
             else
                int_type=>ci%int_type_mv
             endif
             if(argv(3:)=='NN') then
                int_type=MV_INT_TYPE_NN
             elseif(argv(3:)=='LIN') then
                int_type=MV_INT_TYPE_L2P
             elseif(argv(3:)=='SINC') then
                int_type=MV_INT_TYPE_SINC
             elseif(argv(3:)=='PV') then
                int_type=HIST_INT_TYPE_PV
             else
                k=>int_type
             endif
          case('o')
             s=>ci%fout
             cycle
          case('O')
             s=>ci%fout
             mana = ior(mana, CHG_EXTENSION)
             cycle
          case('R')
             x=>ci%red(1)
             ci%red(2)=-1
             ci%red(3)=-1
             mana=ior(mana, REDUCE_R)
          case('b')
             mana = ior(mana, USE_ADAPBIN_FILTER_XY)
             if(len==2) then
                ci%n_abf(1:2)=1
             else
                k=>ci%n_abf(1)
             end if
          case('B')   
             mana = ior(mana, USE_ADAPBIN_FILTER)
             if(len==2) then
                ci%n_abf=1
             else
                k=>ci%n_abf(3)
             endif
          case('s')
             mana = ior(mana, USE_MEDIAN_FILTER_XY)
             if(len==2) then
                ci%n_mf(1:2)=1
             else
                k=>ci%n_mf(1)
             endif
          case('S')
             mana = ior(mana, USE_MEDIAN_FILTER);;; 
             if(len==2) then
                ci%n_mf=1
             else
                k=>ci%n_mf(3)
             endif
          case default
             do j=2,len
                select case(argv(j:j))
                case('W')
                   mana = ior(mana, COUPLED)
                case('a')
                   mana = ior(mana, ELASTIC)
                case('+')
                   mana = ior(mana, IGNORE_NEG)
                case('E')
                   mana = ior(mana, EXTRACT_CCS)
                case('f')
                   mana = ior(mana, IGNORE_CF)
                case('t')
                   mana = ior(mana, SAVE_MOVED)
                case('T')
                   mana = ior(mana, ior(SAVE_MOVED, MATCH_DIMENSION))
                case('v')
                   mana = ior(mana, VERBOSE)
                case('h')
                   call print_help
                   call exit(0)
                case('c')
                   mana = ior(mana, EXCHANGE_INPUT)
                case('Q')
                   mana = ior(mana, NO_OUTPUT)
                case('w')
                   mana = ior(mana, DO_EXP)
                case('e')
                   mana = ior(mana, USE_HISTEQ)
                case('i')
                   mana = ior(mana, GET_MI)
                case('U')
                   mana = ior(mana, MULTIPLY)
                case('V')
                   call print_version
                   call exit(0)
                case default
                   call echo_red("*** Unrecognized option: '"//&
                        trim(argv(2:))//"'")
                end select
             end do
          end select
       else
          ! implicit rule of reading non-optionals 
          if(ci%fin=="") then
             s=>ci%fin
          elseif(ci%fref=="") then
             s=>ci%fref   
             mana = ior(mana, REF_READ)
          elseif(ci%fout=="") then
             s=>ci%fout
          else
             call echo_red("*** Too many arguments at: '"//trim(argv)//"'") 
             return
          endif
       endif
       do_cycle=read_arg()
    end do
    if(cd/=0) call echo_blue("lift off.")
    call echo_blue(">>> This is mana to ground control, mc is "//&
         argv(1:d2hexstr(mana,argv))//&
         ", green?")
    contains 
      logical*4 function read_arg()
        integer*4 ii
        read_arg=.true.
        if(associated(x)) then
           if(len>=3) read(argv(3:len),*,iostat=ii) x
           nullify(x)
        else if(associated(k)) then
           if(len>=3) read(argv(3:len),*,iostat=ii) k
           nullify(k)
        else if(associated(s)) then
           s=argv(1:len)
           ii=0
           nullify(s)
        else 
           read_arg=.false.
           return
        endif
        if(ii/=0) &        
             call echo_red("*** Unrecognized option: '"//trim(argv(2:))//"'")
      end function read_arg
  end subroutine me_first_mana
  !
  subroutine gc_veto(ci)
    type(cmd_info),intent(inout)::ci
    character(PATH_LEN_MAX) str
    if(ci%fin=="") call echo_red("Input file missing.")
    if(green(mana, EXCHANGE_INPUT)) then
       str=ci%fin
       ci%fin=ci%fref
       ci%fref=str
    endif
    select case(ci%ff)
    case(FF_ANALYZE,FF_SHIMADZU)
    case default
       call echo_red("Invalid file format selected.")
    end select
    select case(ci%int_kind)
    case(1,2,4,8,32,64)
    case default
       call echo_red("Invalid internal data type.")
    end select
    call echo_blue("+++ Ground control veto, green.")
  end subroutine gc_veto
  !
  subroutine echo_green(ci)
    type(cmd_info),intent(inout)::ci
    if(ci%red(2)==-1) ci%red(2)=ci%red(1)
    if(ci%red(3)==-1) ci%red(3)=ci%red(1)
    if(ci%d(1)/=0.and.ci%d(2)==0.and.ci%d(3)==0) ci%d(2:3)=ci%d(1)
    if(ci%n_mf(3)/=0) then
       ci%n_mf(1:2)=ci%n_mf(3)
    elseif(ci%n_mf(1)/=0)then
       ci%n_mf(2)=ci%n_mf(1)
    end if
    if(ci%n_abf(3)/=0) then
       ci%n_abf(1:2)=ci%n_abf(3)
    elseif(ci%n_abf(1)/=0)then
       ci%n_abf(2)=ci%n_abf(1)
    end if
    ci%dt=ci%dt/180*PI
    if(ci%fout=="") then
       if(ci%ff==0) then
          ci%fout=trim(ci%fin)//".out"
       else
          ci%fout=trim(ci%fin)
       end if
    else if(green(mana, CHG_EXTENSION)) then
       ci%fout=trim(ci%fin)//"."//trim(ci%fout)
    endif
    call echo_blue("<<< This is ground control to uss long shot, green.")
  end subroutine echo_green
  !
  subroutine echo_blue(str)
    character*(*),intent(in)::str
    if(iand(mana, VERBOSE)/=0) call messn(str)
  end subroutine echo_blue
  !
  subroutine echo_red(str,code)
    use cuty
    use ana
    character*(*),intent(in)::str
    integer*4,intent(in),optional::code
    character*(*),parameter::swan_song="*** mission aborted."
    character*128 buf
    integer*4 l
    if(present(code)) then
       call mess(str)
       l=conc(buf,conc(buf,1,&
            & "code = "),code,1)
       call messn(buf(1:l-1)) 
       if(iand(mana, VERBOSE)/=0) call err_ana(code)
       if(code/=0)then
          call messn(swan_song)
          call exit(code)
       endif
    else
       call messn(str)
    endif
    stop swan_song
  end subroutine echo_red
  !
  logical*4 function green(code,obj)
    integer*4,intent(in)::code
    integer*4,intent(in)::obj
    if(iand(code,obj)==0) then
       green=.false.
    else
       green=.true.
    endif
  end function green
  !
  integer*4 function res_adf(code)
    integer*4,intent(in)::code
    res_adf=ishft(iand(code,ADF_MASK),-24)
  end function res_adf
end module Ground_Control
!=============================================================================!
!=============================================================================!
module uss_long_shot
  use Ground_Control
  implicit none
  real*8,parameter,private::errrel_min=1.d-6 ! <<< 
contains
  recursive subroutine mana_lab(icode,ci,in,out,ref)
    integer*4,intent(in)::icode
    type(cmd_info),intent(in)::ci
    type(ana_data),intent(out)::in
    type(ana_data),intent(out),optional::out
    type(ana_data),intent(in),optional::ref
    type(ana_data) ad
    type(rigid_transform_param) rtp
    type(Interpolation_Data) ipd
    type(gmwr_Data) gmwrd
    integer*4,allocatable::img_i4(:,:,:)
    real*8,allocatable::h2(:,:,:)
    real*8,allocatable::h(:)
    real*8,allocatable::kern(:,:,:)
    real*8 s(3),ss
    integer*4 n(3)
    character*256 str
    integer*4 l,i
    integer*4 rc

    select case(icode)
    case(MANA_MUL_IMG)
       call echo_blue("Multiplying image...")
       select case(in%id%datatype)
       case(DT_IKIND_1)
          ss=sum(in%id%img_i1)
          in%id%img_i1=in%id%img_i1*out%id%img_i1
          ss=ss/sum(in%id%img_i1)
          in%id%img_i1=in%id%img_i1*ss
       case(DT_IKIND_2)
          ss=sum(in%id%img_i2)
          in%id%img_i2=in%id%img_i2*out%id%img_i2
          ss=ss/sum(in%id%img_i2)
          in%id%img_i2=in%id%img_i2*ss
       case(DT_IKIND_4)
          ss=sum(in%id%img_i4)
          in%id%img_i4=in%id%img_i4*out%id%img_i4
          ss=ss/sum(in%id%img_i4)
          in%id%img_i4=in%id%img_i4*ss
       case(DT_RKIND_4)
          ss=sum(in%id%img_r4)
          in%id%img_r4=in%id%img_r4*out%id%img_r4
          ss=ss/sum(in%id%img_r4)
          in%id%img_r4=in%id%img_r4*ss
       case(DT_RKIND_8)
          ss=sum(in%id%img_r8)
          in%id%img_r8=in%id%img_r8*out%id%img_r8
          ss=ss/sum(in%id%img_r8)
          in%id%img_r8=in%id%img_r8*ss
       end select

    case(MASK)
       call echo_blue("Masking...")
       select case(in%id%datatype)
       case(DT_IKIND_1)
          call mask_img(in%id%plo,in%id%pup,in%id%img_i1,&
               &     PVAL_OUT_OF_SIGHT,PVAL_OUT_OF_SIGHT)
       case(DT_IKIND_2)
          call mask_img(in%id%plo,in%id%pup,in%id%img_i2,&
               &     PVAL_OUT_OF_SIGHT,PVAL_OUT_OF_SIGHT)
       case(DT_IKIND_4)
          call mask_img(in%id%plo,in%id%pup,in%id%img_i4,&
               &     PVAL_OUT_OF_SIGHT,PVAL_OUT_OF_SIGHT)
       case(DT_RKIND_4)
          call mask_img(in%id%plo,in%id%pup,in%id%img_r4,&
               &     PVAL_OUT_OF_SIGHT,PVAL_OUT_OF_SIGHT)
       case(DT_RKIND_8)
          call mask_img(in%id%plo,in%id%pup,in%id%img_r8,&
               &     PVAL_OUT_OF_SIGHT,PVAL_OUT_OF_SIGHT)
       end select
       
    case(USE_MEDIAN_FILTER,USE_MEDIAN_FILTER_XY)
       call echo_blue("Calling median filter...")
       select case(in%id%datatype)
       case(DT_IKIND_1)
          call median_filter(in%id%n,in%id%img_i1,ci%n_mf)
       case(DT_IKIND_2)
          call median_filter(in%id%n,in%id%img_i2,ci%n_mf)
       case(DT_IKIND_4)
          call median_filter(in%id%n,in%id%img_i4,ci%n_mf)
       case(DT_RKIND_4)
          call median_filter(in%id%n,in%id%img_r4,ci%n_mf)
       case(DT_RKIND_8)
          call median_filter(in%id%n,in%id%img_r8,ci%n_mf)
       end select
       
    case(USE_ADAPBIN_FILTER,USE_ADAPBIN_FILTER_XY)
       call echo_blue("Calling adapbin filter...")
       select case(in%id%datatype)
       case(DT_IKIND_1)
          call adapbin_filter(in%id%n,in%id%img_i1,ci%n_abf)
       case(DT_IKIND_2)
          call adapbin_filter(in%id%n,in%id%img_i2,ci%n_abf)
       case(DT_IKIND_4)
          call adapbin_filter(in%id%n,in%id%img_i4,ci%n_abf)
       case(DT_RKIND_4)
          call adapbin_filter(in%id%n,in%id%img_r4,ci%n_abf)
       case(DT_RKIND_8)
          call adapbin_filter(in%id%n,in%id%img_r8,ci%n_abf)
       end select
       
    case(USE_HISTEQ)
       call echo_blue("Calling eqhist_img...")
       select case(in%id%datatype)
       case(DT_IKIND_1)
          call eqhist_img(in%id%n,in%id%img_i1,&
               &in%id%plo,in%id%pup,in%id%plo,in%id%pup,in%id%ns)
       case(DT_IKIND_2)
          call eqhist_img(in%id%n,in%id%img_i2,&
               &in%id%plo,in%id%pup,in%id%plo,in%id%pup,in%id%ns)
       case(DT_IKIND_4)
          call eqhist_img(in%id%n,in%id%img_i4,&
               &in%id%plo,in%id%pup,in%id%plo,in%id%pup,in%id%ns)
       case(DT_RKIND_4)
          call eqhist_img(in%id%n,in%id%img_r4,&
               &in%id%plo,in%id%pup,in%id%plo,in%id%pup,in%id%ns)
       case(DT_RKIND_8)
          call eqhist_img(in%id%n,in%id%img_r8,&
               &in%id%plo,in%id%pup,in%id%plo,in%id%pup,in%id%ns)
       end select
       
    case(USE_REMAP)
       call echo_blue("Calling remaphist_img...")
       select case(in%id%datatype)
       case(DT_IKIND_1)
          call remaphist_img(in%id%img_i1,in%id%plo,in%id%pup,in%id%plo,&
               dble(ci%ns))
       case(DT_IKIND_2)
          call remaphist_img(in%id%img_i2,in%id%plo,in%id%pup,in%id%plo,&
               dble(ci%ns))
       case(DT_IKIND_4)
          call remaphist_img(in%id%img_i4,in%id%plo,in%id%pup,in%id%plo,&
               dble(ci%ns))
       case(DT_RKIND_4)
          call remaphist_img(in%id%img_r4,in%id%plo,in%id%pup,in%id%plo,&
               dble(ci%ns))
       case(DT_RKIND_8)
          call remaphist_img(in%id%img_r8,in%id%plo,in%id%pup,in%id%plo,&
               dble(ci%ns))
       end select
       in%c=dble(ci%ns)/in%id%pup
       
    case(CLUSTERIZE)
       call echo_blue("Clusterizing...")
       allocate(img_i4(in%id%n(1),in%id%n(2),in%id%n(3)))
       select case(in%id%datatype)
       case(DT_IKIND_1)
          call kmclusterize_img(ci%nc,in%id%pmin,in%id%pmax,&
               errrel_min,in%id%n,in%id%img_i1,img_i4)
       case(DT_IKIND_2)
          call kmclusterize_img(ci%nc,in%id%pmin,in%id%pmax,&
               errrel_min,in%id%n,in%id%img_i2,img_i4)
       case(DT_IKIND_4)
          call kmclusterize_img(ci%nc,in%id%pmin,in%id%pmax,&
               errrel_min,in%id%n,in%id%img_i4,img_i4)
       case(DT_RKIND_4)
          call kmclusterize_img(ci%nc,in%id%pmin,in%id%pmax,&
               errrel_min,in%id%n,in%id%img_r4,img_i4)
       case(DT_RKIND_8)
          call kmclusterize_img(ci%nc,in%id%pmin,in%id%pmax,&
               errrel_min,in%id%n,in%id%img_r8,img_i4)
       end select
       call echo_blue("Extracting segment...")
       call init_ana_data(ad)
       call cp_ana_data(in,ad)
       call alloc_image_data(ad%id,in%id%n,in%id%datatype)
       do i=1,ci%nc
          select case(in%id%datatype)
          case(DT_IKIND_1)
             where(img_i4==i)
                ad%id%img_i1=in%id%img_i1
             elsewhere
                ad%id%img_i1=in%id%plo
             end where
          case(DT_IKIND_2)
             where(img_i4==i)
                ad%id%img_i2=in%id%img_i2
             elsewhere
                ad%id%img_i2=in%id%plo
             end where
          case(DT_IKIND_4)
             where(img_i4==i)
                ad%id%img_i4=in%id%img_i4
             elsewhere
                ad%id%img_i4=in%id%plo
             end where
          case(DT_RKIND_4)
             where(img_i4==i)
                ad%id%img_r4=in%id%img_r4
             elsewhere
                ad%id%img_r4=in%id%plo
             end where
          case(DT_RKIND_8)
             where(img_i4==i)
                ad%id%img_r8=in%id%img_r8
             elsewhere
                ad%id%img_r8=in%id%plo
             end where
          end select
          l=d2str(i,str)
          ad%fname=trim(ci%fout)//".seg"//str(1:l)
          call cp_image_data(in%id,ad%id)
          call set_image_data(ad%id)
          call print_ana_data(ad)
          call write_ana(ad,ci%out_kind,rc)
       enddo
       select case(in%id%datatype)
       case(DT_IKIND_1)
          in%id%img_i1=img_i4
       case(DT_IKIND_2)
          in%id%img_i2=img_i4
       case(DT_IKIND_4)
          in%id%img_i4=img_i4
       case(DT_RKIND_4)
          in%id%img_r4=img_i4
       case(DT_RKIND_8)
          in%id%img_r8=img_i4
       end select
       deallocate(img_i4)
       ! here we must discard lo and up
       in%id%plo=-HUGE(0.d0)
       in%id%pup=HUGE(0.d0)
       call set_image_data(in%id)
       call print_image_data(in%id,"Segmented")
       call dealloc_image_data(ad%id)
       
    case(EXTRACT_CCS)
       call echo_blue("Extracting connected componets...")
       call init_ana_data(out)
       call cp_ana_data(in,out)
       call label2d_img(in%id,out%id,i)
       l=conc(str,conc(str,1,&
            & "kmax = "),i)
       call messn(str(1:l-1))
       call sort_label2d_img(out%id%n,out%id%img_i4)
       
    case(REDUCE,REDUCE_R,REDUCE_V,REDUCE_AO)
       call echo_blue("Reducing...")
       if(icode/=REDUCE_AO) then
          call init_ana_data(out)
          call cp_ana_data(in,out)
          if(icode==REDUCE_V) then
             out%id%d=max(in%id%d,max(ci%d,0.d0))
          else
             out%id%d=ci%red*in%id%d
          endif
       endif
       call rebin_img(in%id,out%id)
       if(out%id%n(1)==0) call echo_red("*** Memory allocation error.")

    case(GET_IP)
       call echo_blue("Calling infopot filter...")
       call init_ana_data(out)
       call cp_ana_data(in,out)
       call cp_image_data(in%id,out%id)
       call alloc_image_data(out%id,out%id%n,out%id%datatype)
       select case(in%id%datatype)
       case(DT_IKIND_1)
          call infopot_img(5.d0,idnint(in%id%pup-in%id%plo+0.5d0),&
               &       in%id%n,in%id%img_i1,out%id%img_i1)
       case(DT_IKIND_2)
          call infopot_img(5.d0,idnint(in%id%pup-in%id%plo+0.5d0),&
               &       in%id%n,in%id%img_i2,out%id%img_i2)
       case(DT_IKIND_4)
          call infopot_img(5.d0,idnint(in%id%pup-in%id%plo+0.5d0),&
               &       in%id%n,in%id%img_i4,out%id%img_i4)
       case(DT_RKIND_4)
          call infopot_img(5.d0,idnint(in%id%pup-in%id%plo+0.5d0),&
               &       in%id%n,in%id%img_r4,out%id%img_r4)
       case(DT_RKIND_8)
          call infopot_img(5.d0,idnint(in%id%pup-in%id%plo+0.5d0),&
               &       in%id%n,in%id%img_r8,out%id%img_r8)
       end select
       
    case(MOVE)
       call echo_blue("Moving image...")
       call init_ana_data(out)
       call cp_ana_data(in,out)
       if(green(mana,MATCH_DIMENSION).and.present(ref)) then          
          call cp_image_data(ref%id,out%id)
       else
          call cp_image_data(in%id,out%id)
       end if
       call init_ipd(ipd,ci%int_type_mv,&
            nx_lnpt=ci%int_type_mv,&
            nx_sinct=3,d_sinct=128)
       rtp%d=ci%dx
       rtp%t=ci%dt
       rtp%m=ci%dm
       rtp%s=ci%ds
       call mv_img(rtp,ci%int_type_mv,in%id,out%id)
       call uninit_ipd(ipd)
       
    case(GET_MI)
       call echo_blue("Evaluating mi...")
       rtp%d=0
       rtp%t=0
       if(ci%nh/=0)then
          in%id%ns=ci%nh
          out%id%ns=ci%nh
          call set_image_data(in%id,flag=SIDF_SET_NSAB)
          call set_image_data(out%id,flag=SIDF_SET_NSAB)
       endif
       if(allocated(h2)) deallocate(h2)
       allocate(h2(in%id%ns,out%id%ns,2))
       call init_gmwr(gmwrd,ci%int_type_in,&
            nx_lnpt=ci%int_type_in,&
            nx_sinct=3,d_sinct=128,pdf_smtw=ci%pdf_smtw)
       call regist_buty_work_space(h2)

       in%id%H=get_H(in%id%n,in%id%img_r8,in%id%ns,in%id%a,in%id%b)
       out%id%H=get_H(out%id%n,out%id%img_r8,out%id%ns,out%id%a,out%id%b)

       l=conc(str,conc(str,1,&
            & "mi = "),get_mi_w_rtp(rtp,in%id,out%id))
       call messn(str(1:l-1))
       call uninit_gmwr(gmwrd)
       deallocate(h2)
       return
       
    case(CUT_BG)
       call echo_blue("Searching histogram peak...")
       allocate(h(in%id%ns))
       call get_hist(in%id%n,in%id%img_r8,in%id%plo,in%id%pup,in%id%ns,h)
       i=find_bg_level(in%id%ns,h,ci%bg_level)
       deallocate(h)
       if(i<=0) call echo_red("find_bg_level failed.")
       in%id%plo=(i-in%id%b)/in%id%a
       call mana_lab(MASK,ci,in)
       return
       
    case(USE_GSMT)
       call echo_blue("Calling Gaussian smoothing...")
       s=ci%gsmt_s/in%id%d
       call get_gmask_size(s,1.d0,n)
       allocate(kern(-n(1):n(1),-n(2):n(2),-n(3):n(3)))
       call get_gmask(s,n,kern)
       select case(in%id%datatype)
       case(DT_IKIND_1)
          call filter_img(in%id%n,n,kern,in%id%img_i1)
       case(DT_IKIND_2)
          call filter_img(in%id%n,n,kern,in%id%img_i2)
       case(DT_IKIND_4)
          call filter_img(in%id%n,n,kern,in%id%img_i4)
       case(DT_RKIND_4)
          call filter_img(in%id%n,n,kern,in%id%img_r4)
       case(DT_RKIND_8)
          call filter_img(in%id%n,n,kern,in%id%img_r8)
       end select
    end select
    
    if(present(out)) then
       call set_image_data(out%id)
       if(green(mana, VERBOSE)) call print_image_data(out%id)
    else
       call set_image_data(in%id)
       if(green(mana, VERBOSE)) call print_image_data(in%id)
    end if

  end subroutine mana_lab
end module uss_long_shot
!=============================================================================!
program faely
!
! Projects BEGIN: 2004/04/05
! Tools for medical images in analyze format.
!
! Current mission:
!   Construnction of a model for arbitrary image-to-image matching and find 
!   its unique solution.
! 
! Past missions: (newer to older)
!   1) image segmentation => connected component analysis (2D) -> done. 
!   2) rigid body transformation model for image matching. -> done.
!   3) multiresolution handling (rebinnig in 3D). -> done.
!   4) n-point Laglange interpolation (3D) -> done.
!   5) completely free deformation of a image -> aborted.
!
! Bibliography:
!   1) "A fully automatic multimodality image registration algorithm",
!      Babak A. et. al., Journal of Computer Assisted Tomography, 19(1995)615.
!   2) "PET-CT image registration in the chest using free-form deformations",
!      Mattes D. et. al., IEEE trans. med. img. 22(2003)120.
!   3) "An overlap invariant entropy measure of 3D medical image alignment",
!      C. Studholme, et. al., Pattern Recogn. 32(1999)71.
!   4) "Image matching as a diffusion process: an analogy with Maxwell's 
!       demons", J.-P. Thirion, Medical Image Analysis 2(1998)243.
!   5) "Information theoretic clustering", Erhan Gokcay et. al., IEEE trans. 
!      Pattern Analysis and Machine Intelligence 20(2002)158.
!
  use Ground_Control
  use uss_long_shot
  use puty
  implicit none
  type(cmd_info) ci
  type(ana_data),pointer::ref,tar
  type(ana_data),pointer::rref,ttar 
  type(ana_data),pointer::reff,tarr
  type(ana_data) orig
  type(rigid_transform_param),target::rtp
  type(gmwr_data) gmwrd
  real*8 errabs_min,errrelobj_min,errrel_min
  real*8 errabs,errrel,errrelobj
  real*8 lambda
  integer*4 rc,n
  real*8,target::x(12),dx(12),xmin(12),xmax(12),xout(12)
  real*8,pointer::x_d(:),x_r(:),dx_r(:),dx_d(:),xout_d(:),xout_r(:)
  real*8,pointer::xmin_d(:),xmax_d(:),xmin_r(:),xmax_r(:)
  real*8,pointer::x_s(:),x_m(:),dx_s(:),dx_m(:),xout_s(:),xout_m(:)
  real*8,pointer::xmin_s(:),xmax_s(:),xmin_m(:),xmax_m(:)
  real*8 afn(3,4)
  character*256 str
  real*8 time0,time1
  real*8 dmin,mi
  !
  call me_first_mana(ci)
  call gc_veto(ci)
  call echo_blue("Reading image: "//trim(ci%fin))
  allocate(tar)
  rc=res_adf(mana)
  select case(ci%ff)
  case(FF_ANALYZE)
     call read_ana(ci%fin,tar,ci%int_kind,rc)
     call read_ana(ci%fin,orig,ci%int_kind,rc)
  case(FF_SHIMADZU)
     call read_shimadzu(ci%fin,tar,rc)
  end select
  if(rc/=0) call echo_red("Read_ana retruned with ",rc)
  tar%id%pup=ci%upper
  tar%id%plo=ci%lower
  !
  ORIG%id%pup=ci%upper
  ORIG%id%plo=ci%lower
  !
  call set_image_data(tar%id)
  call set_image_data(ORIG%id)
  call print_ana_data(ORIG,"Original")
  call print_ana_data(tar,"Target")
  ttar=>tar
  !
  if(green(mana, REF_READ)) then
     allocate(ref)
     call echo_blue("Reading image: "//trim(ci%fref))
     rc=res_adf(mana)
     call read_ana(ci%fref,ref,ci%int_kind,rc)
     if(rc/=0) call echo_red("Read_ana retruned with ",rc)
     ref%id%pup=ci%upper
     ref%id%plo=ci%lower
     call set_image_data(ref%id)
     call print_ana_data(ref,"Reference")
     rref=>ref
  else
     rref=>tar
  endif
  !
  ! Mission start
  call echo_green(ci)
  if(rref%id%n(3)==ttar%id%n(3)) then
     rref%id%d(3)=min(rref%id%d(3),ttar%id%d(3))
     ttar%id%d(3)=rref%id%d(3)
  end if
  !
  if(green(mana, MOVE)) then
     allocate(rref%ad)
     call mana_lab(MOVE,ci,rref,rref%ad)
     if(green(mana, REF_READ)) then
        rref=>rref%ad
     else
        ttar=>rref%ad
     end if
  end if
  !
  if(green(mana, USE_GSMT)) then
     call mana_lab(USE_GSMT,ci,ttar)
     if(green(mana, REF_READ)) call mana_lab(USE_GSMT,ci,rref)
  endif
  !
  if(green(mana, CUT_BG)) then
     call mana_lab(CUT_BG,ci,ttar)
     if(green(mana, REF_READ)) call mana_lab(CUT_BG,ci,rref)
  end if
  !
  if(green(mana, MASK)) then
    call mana_lab(MASK,ci,ttar)
    if(green(mana, REF_READ)) call mana_lab(MASK,ci,rref)
  endif
  !  
  if(green(mana, ior(USE_MEDIAN_FILTER,USE_MEDIAN_FILTER_XY))) then
     call mana_lab(USE_MEDIAN_FILTER,ci,ttar)
     if(green(mana, REF_READ)) call mana_lab(USE_MEDIAN_FILTER,ci,rref)
  endif
  !
  if(green(mana, USE_HISTEQ)) then
     call mana_lab(USE_HISTEQ,ci,ttar)
     if(green(mana, REF_READ)) call mana_lab(USE_HISTEQ,ci,rref)
  endif
  !
  if(green(mana, USE_REMAP)) then
     call mana_lab(USE_REMAP,ci,ttar)
     if(green(mana, REF_READ)) call mana_lab(USE_REMAP,ci,rref)
  endif
  !
  if(green(mana, CLUSTERIZE)) then
     call mana_lab(CLUSTERIZE,ci,ttar)
     if(green(mana, REF_READ)) call mana_lab(CLUSTERIZE,ci,rref)
  endif
  !
  if(green(mana, ior(USE_ADAPBIN_FILTER,USE_ADAPBIN_FILTER_XY))) then
     call mana_lab(USE_ADAPBIN_FILTER,ci,ttar)
     if(green(mana, REF_READ)) call mana_lab(USE_ADAPBIN_FILTER,ci,rref)
  endif
  !
  nullify(tarr,reff)
  if(green(mana, REDUCE)) then
     tarr=>ttar
     allocate(ttar%ad)
     call mana_lab(iand(mana,REDUCE),ci,ttar,ttar%ad)
     ttar=>ttar%ad
     if(iand(mana, REF_READ)/=0) then
        reff=>rref
        allocate(rref%ad)
        call mana_lab(iand(mana,REDUCE),ci,rref,rref%ad)
        rref=>rref%ad
     else
        rref=>ttar
     endif
  end if
  !
  if(green(mana, EXTRACT_CCS)) then
     allocate(ttar%ad)
     call mana_lab(EXTRACT_CCS,ci,ttar,ttar%ad)
     ttar=>ttar%ad
  end if
  !
  if(green(mana, GET_IP)) then
     allocate(ttar%ad)
     call mana_lab(GET_IP,ci,ttar,ttar%ad)
     ttar%ad%fname="out2"
     call write_ana(ttar%ad,ci%out_kind,rc)
     ttar=>ttar%ad
  endif
  !
  if(green(mana, GET_MI)) call mana_lab(GET_MI,ci,ttar,rref)
  !
  if(green(mana, MULTIPLY)) call mana_lab(MANA_MUL_IMG,ci,ttar,rref)
  !
  if(.not.green(mana, DO_EXP)) goto 999
  !
  !
  !
  call echo_blue("Matching BEGIN...")
  call set_image_data(rref%id,flag=ior(SIDF_ADJ_LOUP,SIDF_SET_ALL))
  call set_image_data(ttar%id,flag=ior(SIDF_ADJ_LOUP,SIDF_SET_ALL))
  call print_image_data(rref%id,"Reference")
  call print_image_data(ttar%id,"Target")
  !
  call set_pair(ttar%id,rref%id)
  call set_rtp(rtp)
  call init_gmwr(gmwrd,ci%int_type_in,&
       nx_lnpt=ci%int_type_in,&
       nx_sinct=3,d_sinct=128,pdf_smtw=ci%pdf_smtw)
  ! pointer asignments and all
  xmin_d=>xmin(1:3)
  xmin_r=>xmin(4:6)
  xmin_m=>xmin(7:9)
  xmin_s=>xmin(10:12)
  !
  xmax_d=>xmax(1:3)
  xmax_r=>xmax(4:6)
  xmax_m=>xmax(7:9)
  xmax_s=>xmax(10:12)
  !
  x_d=>x(1:3)
  x_r=>x(4:6)
  x_m=>x(7:9)
  x_s=>x(9:12)
  !
  xout_d=>xout(1:3)
  xout_r=>xout(4:6)
  xout_m=>xout(7:9)
  xout_s=>xout(10:12)
  !
  dx_d=>dx(1:3)
  dx_r=>dx(4:6)
  dx_m=>dx(7:9)
  dx_s=>dx(10:12)
  !
  ! init
  x=0
  rtp%m=1
  rtp%s=0
  errrelobj=0
  errabs=0
  !==========================================================================
  ! translation 
  !
  xmax_d=max(rref%id%xmax,ttar%id%xmax)/4
  xmin_d=min(rref%id%xmin,ttar%id%xmin)/4
  xmax_d(3)=ttar%id%xmax(3)-rref%id%xmax(3)/2
  xmin_d(3)=ttar%id%xmin(3)-rref%id%xmin(3)/2
  errabs_min=min(minval(ttar%id%d),minval(rref%id%d))/2
  errrel_min=1.d-4
  errrelobj_min=1.d-4
  rtp%d=x_d
  rtp%t=x_r
  call mess("SIMPLEX D: ")
  call cpu_time(time0)
  call simplex(Cost_D,3,rtp%d,xmin_d-rtp%d,xmax_d-rtp%d,&
       & errabs_min,errrel_min,0.d0,errabs,errrel,errrelobj,xout_d,&
       & num_ite=n,minobj_out=mi)
  call cpu_time(time1)
  rtp%d=xout_d
  call print_results(xout_d)
  !==========================================================================
  ! rotation
  !
  xmin_d=xmin_d/2
  xmax_d=xmax_d/2
  xmin_r=-piq
  xmax_r=+piq
  dx_r=min(sqrt(sum(ttar%id%d**2)/dot_product(ttar%id%d,ttar%id%n)**2),&
       sqrt(sum(rref%id%d**2)/dot_product(rref%id%d,rref%id%n)**2))
  errabs_min=0
  errrel_min=errrel*1.d-2
  errrelobj_min=errrelobj*1.d-2
  x_d=rtp%d
  x_r=rtp%t
  x_r(1)=-16./180.*pi
  if(green(mana, COUPLED)) then
     call mess("SIMPLEX RD: ")
     call cpu_time(time0)
     call simplex(Cost_RD,6,x(1:6),xmin(1:6)-x(1:6),xmax(1:6)-x(1:6),&
          & errabs_min,errrel_min,errrelobj_min,&
          & errabs,errrel,errrelobj,xout(1:6),&
          & num_ite=n,minobj_out=mi)
     rtp%d=xout_d
  else
     call mess("SIMPLEX R: ")
     call cpu_time(time0)
     call simplex(Cost_R,3,x_r,xmin_r-x_r,xmax_r-x_r,&
          & errabs_min,errrel_min,errrelobj_min,&
          & errabs,errrel,errrelobj,xout_r,&
          & num_ite=n,minobj_out=mi)
  endif
  call cpu_time(time1)
  rtp%t=xout_r
  call print_results(xout_d,xout_r/pi*180)
  !==========================================================================
  ! magnification 
  !
  if(green(mana, ELASTIC)) then
     x_m=1
     rtp%s=0
     xmin_m=-0.5d0
     xmax_m= 1.5d0
     errabs_min=0
     errrel_min=max(1.d-4,errrel)       ! of cost_r
     errrelobj_min=max(1.d-4,errrelobj) ! of cost_r
     call mess("SIMPLEX M: ")
     call cpu_time(time0)
     call simplex(Cost_M_nr,3,x_m,xmin_m-x_m,xmax_m-x_m,&
          & errabs_min,errrel_min,errrelobj_min,&
          & errabs,errrel,errrelobj,xout_m,&
          & num_ite=n,minobj_out=mi)
     call cpu_time(time1)
     rtp%m=xout_m
     call print_results(xout_m)
  end if
  !==========================================================================
  ! sheer 
  !
  if(green(mana, ELASTIC)) then
     x_s=0
     xmin_s=-.5d0
     xmax_s= .5d0
     errabs_min=0
     errrel_min=max(1.d-4,errrel)       ! of cost_r
     errrelobj_min=max(1.d-4,errrelobj) ! of cost_r
     call mess("SIMPLEX S: ")
     call cpu_time(time0)
     call simplex(Cost_S_nr,3,x_s,xmin_s-x_s,xmax_s-x_s,&
          & errabs_min,errrel_min,errrelobj_min,&
          & errabs,errrel,errrelobj,xout_s,&
          & num_ite=n,minobj_out=mi)
     call cpu_time(time1)
     rtp%s=xout_s
     call print_results(xout_s)
  end if
  !
  if(associated(tarr).and.associated(reff)) then
     dmin=get_max_d(tarr%id%d,reff%id%d)
     ! change resolution
     if(dmin<get_max_d(ttar%id%d,rref%id%d)) then        
        ttar%id%d=dmin
        rref%id%d=dmin
        call mana_lab(REDUCE_AO,ci,reff,rref)
        call mana_lab(REDUCE_AO,ci,tarr,ttar)
        call set_image_data(rref%id,flag=ior(SIDF_ADJ_LOUP,SIDF_SET_ALL))
        call set_image_data(ttar%id,flag=ior(SIDF_ADJ_LOUP,SIDF_SET_ALL))
        call print_image_data(rref%id,"Reference")
        call print_image_data(ttar%id,"Target")
        call set_pair(ttar%id,rref%id)
     end if
  endif
  !
  dx_d=min(rref%id%d,ttar%id%d)
  lambda=sqrt(sum(dx_d**2))
  errabs_min=0
  errrel_min=0
  errrelobj_min=1.d-5
  call mess("GNEWTON D: ")
  call cpu_time(time0)
  if(green(mana, ELASTIC)) then
     call gnewton(Cost_D_nr,3,rtp%d,dx_d,&
          & lambda,errabs_min,errrel_min,errrelobj_min,&
          & errabs,errrel,errrelobj,xout_d,&
          & num_ite=n,minobj_out=mi)
  else
     call gnewton(Cost_D,3,rtp%d,dx_d,&
          & lambda,errabs_min,errrel_min,errrelobj_min,&
          & errabs,errrel,errrelobj,xout_d,&
          & num_ite=n,minobj_out=mi,rc_out=rc)
  endif
  rtp%d=xout_d
  call cpu_time(time1)
  call print_results(xout_d)
  !
  dx_r=min(sqrt(sum(ttar%id%d**2)/dot_product(ttar%id%d,ttar%id%n)**2),&
       sqrt(sum(rref%id%d**2)/dot_product(rref%id%d,rref%id%n)**2))
  lambda=sqrt(sum(dx_r**2))*10*2
  errabs_min=0
  errrel_min=0
  errrelobj_min=1.d-5
  call mess("GNEWTON R: ")
  call cpu_time(time0)
  if(green(mana, ELASTIC)) then
     call gnewton(Cost_R_nr,3,rtp%t,dx_r,&
          & lambda,errabs_min,errrel_min,errrelobj_min,&
          & errabs,errrel,errrelobj,xout_r,&
          & num_ite=n,minobj_out=mi)
  else
     call gnewton(Cost_R,3,rtp%t,dx_r,&
          & lambda,errabs_min,errrel_min,errrelobj_min,&
          & errabs,errrel,errrelobj,xout_r,&
          & num_ite=n,minobj_out=mi)
  end if
  rtp%t=xout_r
  call cpu_time(time1)
  call print_results(xout_r*180/pi)
  !
  call set_affine(rtp,afn)
  call messn("afn:")
  call print_affine(afn)
  call echo_blue("Optimization END")
  !
  if(green(mana, DUMP_PDF)) then
     call echo_blue("Writing final pdf...")
     call dump_h2d("p12fin.dat",tar%fname,ref%fname)
  endif
  !
  if(green(mana, SAVE_MOVED)) then
     call echo_blue("Writing transformed image...")
     ci%dx=rtp%d
     ci%dt=rtp%t
     ci%dm=rtp%m
     ci%ds=rtp%s
     if(.not.associated(tar%ad)) allocate(tar%ad) ! <<< 
!!$     call dealloc_image_data(tar%id)               ! FIX ME <<<
!!$     call cp_image_data(ORIG%id,tar%id)            ! FIX ME <<<
!!$     call alloc_image_data(tar%id,ORIG%id%n)       ! FIX ME <<<   
!!$     tar%id%img_r8=ORIG%ID%IMG_R8                  ! FIX ME <<<
     call mana_lab(MOVE,ci,ORIG,out=tar%ad,ref=ref)   ! FIX ME <<<
     tar%ad%fname=trim(tar%fname)//".reg"
     call echo_blue("Writing image: "//trim(tar%ad%fname))
     call write_ana(tar%ad,ci%out_kind,rc)
     if(rc/=0) call echo_red("write_ana returned with ",rc)
  endif
999 continue
  if(.not.green(mana, NO_OUTPUT)) then
     call echo_blue("Writing image: "//trim(ci%fout))
     ttar%fname=ci%fout     
     call print_ana_data(ttar)
     call write_ana(ttar,ci%out_kind,rc)
     if(rc/=0) call echo_red("write_ana returned with ",rc)
  endif
contains
  subroutine print_results(results1,results2)
    real*8,intent(in)::results1(3)
    real*8,intent(in),optional::results2(3)
    call messn(str(1:f2str(time1-time0,str))//' sec')
    call messn('obj '//str(1:f2str(mi,str)))
    call messn('n = '//str(1:d2str(n,str)))
    call messn('errabs = '//str(1:f2str(errabs,str)))
    call messn('errrel = '//str(1:f2str(errrel,str)))
    call messn('errrelobj = '//str(1:f2str(errrelobj,str)))
    call messn(str(1:f2str(results1,str)))
    if(present(results2)) call messn(str(1:f2str(results2,str)))
  end subroutine print_results
end program faely
